/* exe_stk_extinction.cpp */

#include <iostream>
#include <string>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <list>

#include <quigon.h>
#include "stk_para.h"
#include "stk_lidar.h"
#include "stk_lidar_io.h"
#include "basic_process.h"
#include "noise_filter.h"
#include "fernald.h"
#include "cloud_detect.h"
#include "output_lidarpro.h"

using namespace std;
using namespace quigon;
using namespace boost::posix_time;

int main(int argc, char *argv[])
{
	string start_t, end_t;
	double distance;
	if (argc < 4)
	{
		cout << "usage: stk_extinction yyyymmddThhmmss yyyymmddThhmmss distance" << endl;
		cout << "The capital 'T' in the 'iso time string' is optional" << endl;
		exit(0);
	}
	else
	{
		start_t = argv[1];
		end_t = argv[2];
		distance = atof(argv[3]);
	}

	// Check the date_time format, if not iso string, add 'T' into it.
	if (start_t[8] != 'T')
		start_t.insert(8, 1, 'T');
	if (end_t[8] != 'T')
		end_t.insert(8, 1, 'T');

	// Load lidar parameters
	STK_Para para("../settings/stk_lidar_para.txt");

	int max_line_size = para.distance_i[distance] + 1;

	// Load proper temperature and pressure "profiles"
	// to calculate backscatter coefficient of atmosphere molecules: beta_m
	Atmos_Profile p_prof = load_atmos_profile("../settings/atmos_profiles/US_1976.p_profile", para);
	Atmos_Profile t_prof = load_atmos_profile("../settings/atmos_profiles/US_1976.t_profile", para);
	Atmos_Profile beta_m = cal_beta_m(t_prof, p_prof, 523.0, 0.0284);

	// Loading the priori data
//	cout << "Loading priori data..." << endl;
//	Time_Serie fe_distance("../settings/lidar/stk_priori_distance.txt");
//	Time_Serie fe_value("../settings/lidar/stk_priori_value.txt");

	// Making a container for the data
	list<STK_Line> container;

	// import the data.
	cout << "Importing data..." << endl;
	import(
	        container,
	        from_iso_string(start_t),		// boost::posix_time::from_iso_string
	        from_iso_string(end_t),		// take care of the "T" in the string
	        para
	);

	cout << "lines num: " << container.size() << endl;
	
	for(list<STK_Line>::iterator i = container.begin();
		i != container.end();
		++i)
	{
		i->resize(max_line_size);
	}

	cout << "Making an averager" << endl;
	Average_Machine<STK_Line> aver(60, stk_aver_attr);
	// Make a functor for correction, and more options are available
	// e.g. Correction<STK_Line> my_corr(para, Status::BACKGROUND | Status::AFTERPULSE | Status::DISTANCE);
	// default is doing every possible correction
	cout << "Making a corrector" << endl;
	Correction<STK_Line, STK_Para> corr1(
		para, 
		Status::BACKGROUND | 
		Status::AFTERPULSE |
		Status::OVERLAP    |
		Status::ENERGY, 
		false);
	Correction<STK_Line, STK_Para> corr2(
		para,
		Status::DISTANCE,
		false);

	// Make a functor for filling the lowest unseen area.
	// here 5 means use the average of lowest 5 known data to fill the unknown area.
	cout << "Making a filler" << endl;
	Fill_Lowest filler(para, 5);

	// Make a filter
	cout << "Making a filter ..." << endl;
	Serie<string, double> filter_para("../settings/tune/stk_pbl_detect_filter.txt");
	Noise_Filter_Proto filter(
	        filter_para["low_index"],
	        filter_para["high_index"],
	        filter_para["low_buddy_num"],
	        filter_para["high_buddy_num"],
	        para.zero_meter_index,
	        max_line_size - 1);
	double trigger_value = filter_para["trigger_value"];
	size_t lowest_trigger_index = filter_para["lowest_trigger_index"];

	// Make a functor for extinction retrieval
	// here 6.0 means only calculate to 6.0 km away, to save time.
	cout << "Making a retriever" << endl;
	Fernald<STK_Line, STK_Para> fernald(para, beta_m, distance, -999.9, -999.9, false);

	// On output
	string output_file_name_base = string("../output/stk_") + start_t + "_" + end_t;
	string::iterator str_iter =
	        remove(output_file_name_base.begin(), output_file_name_base.end(), 'T');
	output_file_name_base.erase(str_iter, output_file_name_base.end());

	string extinction_file_name = output_file_name_base + ".extinctn.txt";
	string ex_diff_file_name = output_file_name_base + ".ex_diff.txt";
	string cbh_file_name = output_file_name_base + ".ex_cbh.txt";
	string pbl_file_name = output_file_name_base + ".ex_pbl.txt";

	Output_Lidarpro printer_extinction(extinction_file_name, para.zero_meter_index);
	Output_Lidarpro printer_diff(ex_diff_file_name, para.zero_meter_index);

	// Average, correction, filter, fill, and retrieval
	aver(container);
	for (list<STK_Line>::iterator i = container.begin();
	     i != container.end();
	     ++i)
	{
//		cout << "######################################" << endl;
//		print_line_brief(*i);
//		cout << "--------------------------------------" << endl;
		corr1(*i);
		filter(*i, trigger_value, lowest_trigger_index);
		corr2(*i);
		filler(*i);
//		print_line_brief(*i);

//		cout << "--------------------------------------" << endl;
		fernald(*i);
//		print_line_brief(*i);
	}

	Time_Serie cbh;
	Time_Serie pbl;
	Cloud_Detect cloud_detector("../settings/tune/stk_ex_cloud_detect.para.txt", 
		Cloud_Detect::SINGLE | Cloud_Detect::REPEAT | Cloud_Detect::DIFF);
	Serie<string, double> pbl_detect_para("../settings/tune/stk_ex_pbl_detect.para.txt");
	double pbl_detect_ratio = pbl_detect_para["pbl_detect_ratio"];
	double pbl_detect_thres = pbl_detect_para["pbl_detect_thres"];
	
	list<STK_Line> diff_container(container);
	for (list<STK_Line>::iterator ex_iter = container.begin(), diff_iter = diff_container.begin();
		ex_iter != container.end() && diff_iter != diff_container.end();
		++ex_iter, ++diff_iter)
	{
		int cbh_index = cloud_detector(*ex_iter);
		if(cbh_index != 0)
			cbh.add_entry(ex_iter->t(), para.distance[cbh_index] * (*ex_iter)["sin_elev_angle"]);		
		adjacent_difference(ex_iter->begin(), ex_iter->end(), diff_iter->begin());

		double current_min = 0;
		size_t current_index = 0;
		for(size_t i = 0; i < (cbh_index == 0 ? diff_iter->size() : cbh_index); ++i)
		{
			(*diff_iter)[i] *= exp(-(*ex_iter)[i] * pbl_detect_ratio);
			if((*diff_iter)[i] < current_min)
			{
				current_index = i;
				current_min = (*diff_iter)[i];
			}
		}	
		if(cbh_index == 0 ||
			current_min < pbl_detect_thres)
		{
			pbl.add_entry(diff_iter->t(), para.distance[current_index] * (*diff_iter)["sin_elev_angle"]);
		}
	}


	// fill time gaps.
	cout << "Filling time gap..." << endl;
	fill_time_gap(container, 60);
	fill_time_gap(diff_container, 60);

	cout << "Outputing to file: " << output_file_name_base + " .extinctn.txt/.ex_diff.txt/.ex_cbh.txt " << endl;
	for (list<STK_Line>::iterator i = container.begin();
	     i != container.end();
	     ++i)
	{
		printer_extinction(*i);
	}
	for (list<STK_Line>::iterator j = diff_container.begin();
		j != diff_container.end();
		++j)
	{
		printer_diff(*j);
	}
	cbh.save(cbh_file_name);
	pbl.save(pbl_file_name);
	
	return 0;
}

